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Abstract 

A  three-dimensional  (3D)  parabolic  equation  acoustical  propagation  code  has  been  developed 
and  run  successfully.  The  code  is  written  in  the  MATLAB  language  and  runs  in  the  MATLAB 
environment.  The  code  has  been  implemented  in  two  versions,  applied  to 

(1)  Horizontal  low-frequency  (100  to  500  Hz)  propagation  through  the  shallow 
water  waveguide  environment; 

(2)  Vertical  high-frequency  propagation  (6  to  15  kHz)  to  study  normal-incidence 
reflection  from  the  lower  side  of  the  ocean  surface. 

The  first  edition  of  the  code  reported  on  here  does  not  implement  refinements  that  are  often 
found  in  2D  propagation  models,  such  as  allowing  density  to  vary,  optimally  smoothing  sound- 
speed  discontinuities  at  the  water/seabed  interface,  and  allowing  an  omni-directional  source.  The 
code  is  part  of  a  development  effort  to  test  the  applicability  of  2D  (and  N  by  2D)  models,  which 
have  more  refinements  than  this  model,  to  the  study  of  fully  3D  propagation  problems,  such  as 
sound  transiting  steep  nonlinear  coastal-area  internal  waves  and/or  sloping  terrain,  and  to  provide 
a  numerical  tool  when  the  full  3D  solution  is  needed. 

Introduction 

A  computational  code  has  been  implemented  to  study  truly  three-dimensional  acoustical 
propagation.  The  code  and  some  initial  results  obtained  with  the  code  are  described  here.  The 
method  that  is  employed  is  the  split-step  Fourier  (SSF)  solution  of  the  parabolic  acoustic  wave 
equation  (PE)  [1,2].  The  equation  describes  one-way  solutions  for  waves  propagating  from  a 
monochromatic  source.  So-called  marching  algorithms  can  be  used  to  find  the  steady-state  one¬ 
way  solutions.  The  SSF  technique  divides  propagation  over  each  distance  increment  through  a 
heterogeneous  sound  speed  (refractive  index)  environment  into  step-by-step  “free  space” 
propagation  through  a  medium  having  a  fixed  wavenumber  (usually  related  to  a  fixed  sound 
speed,  often  the  mean),  and  periodically  introduced  (at  each  step)  phase  fluctuations  consistent 
with  departures  from  that  fixed  speed.  The  free  space  propagation  is  handled  in  the  wavenumber 
domain,  and  the  phase  anomalies  are  introduced  in  the  spatial  domain.  Amplitude  effects  such  as 
absorption  are  introduced  with  the  phase  anomalies.  Thus,  each  step  requires  a  Fourier  transform 
and  an  inverse  Fourier  transform,  with  the  spatially  marched  complex  acoustic  field  solution  T 
given  by 

W(x+  S)  =  r'[G-(F[PW(x)])] 

where  P=Ar  exp(  - ik()US )  is  the  phase  (and  amplitude)  operator  in  the  spatial  domain,  G  is  the 
propagator  in  the  wavenumber  domain,  F  is  the  Fourier  transform  operator,  U  =  (c-co )  /  Co ,  Co 


is  a  reference  sound  speed  and  c  is  sound  speed,  ko=  co/co  is  wavenumber,  co  is  frequency,  and  5 
is  the  distance  increment  in  the  x  direction. 

In  ocean  acoustics,  the  PE  is  usually  implemented  in  vertical  slices  (vertical  planes).  At 
each  step,  the  vertical  dimension  and  the  vertical  wavenumber  form  the  transform  variable  pair. 

A  few  2D  SSF  PE  implementations  are  available,  including  the  Monterey-Miami  PE  and  the 
University  of  Miami  PE  which  preceded  it.  These  methods  utilize  single-dimension  fast  Fourier 
transforms  (FFTs).  For  horizontal  propagation,  these  methods  implement  the  surface  boundary 
condition  with  an  out-of  phase  image  source  in  an  image  domain  on  the  opposite  side  of  the 
surface  from  the  water  and  seafloor,  meaning  that  the  numerical  domain  size  is  twice  that  of  the 
modeled  water-column  and  the  modeled  portion  of  the  sea  floor. 

The  extension  of  the  2D  SSF  PE  method  to  a  Cartesian  coordinate  3D  environment 
primarily  involves  changing  from  a  1 D  to  a  2D  Fourier  transform.  This  adds  substantially  to  the 
calculation  time  required  to  cover  a  specific  distance  and  depth  of  water  at  a  specific  grid 
resolution.  In  addition  to  the  added  run  time,  larger  amounts  of  memory  are  also  required  to 
perform  the  calculation.  The  code  reported  on  here  uses  the  aforementioned  image-domain  2-D 
FFT  technique.  The  wide-angle  variant  of  the  propagation  operator  is  used  [2].  The  aspect  ratio 
of  the  vertical  and  horizontal  grid  intervals  is  adjustable,  as  is  the  domain  size,  the  dimensions  of 
the  FFT,  and  the  propagation  step  length.  This  code  originates  from  3D  propagation  codes  used 
for  optical  studies  [3,4],  A  trial  MATLAB  “  *  version  was  written  by  John  Colosi,  from  which 
this  code  was  derived. 

Note  that  a  different  3D  PE  code,  FOR3D,  has  been  used  for  some  years  in  ocean  acoustics, 
perhaps  most  recently  by  Naval  Research  Laboratory  investigators  [5,6],  FOR3D  uses  a  range- 
azimuth  computational  grid  and  implements  narrow-angle  azimuthal  coupling. 

There  are  number  of  issues  still  to  be  resolved  in  order  to  apply  the  technique  in  a  reliable 
fashion  to  ocean  acoustics  problems.  One  unknown  is  the  grid  resolution  needed  to  ensure  a 
good  solution.  The  definition  of  a  source  field  and  source  level  has  not  yet  been  resolved. 

Finally,  the  required  absorption  of  sound  at  the  lateral  edges  is  being  investigated.  Progress  on 
these  in  two  distinct  implementations  is  described  here,  and  example  solutions  are  shown.  The 
two  implementations  are  100-500  Hz  propagation  up  to  10’s  of  kilometers  in  a  shallow  water 
waveguide,  and  o(10  kHz)  upward  then  downward  vertical  propagation  with  a  reflection  at  the 
rough  (wavy)  sea  surface  intended  to  model  signals  from  inverted  echo  sounders.. 

Shallow  Water  Waveguide 

Code  for  modeling  the  coastal  shallow-water  waveguide  has  been  completed.  Some 
preliminary  results  are  shown  here.  All  results  shown  in  this  section  are  for  400-Hz  propagation. 
Grid  sizes  of  1024  x  1024  and  2048  x  2048  have  been  used.  The  propagation  step  used  here  is 
two  acoustic  wavelengths  (X)  in  length  (7.5  m  in  this  case).  The  terminology  here  is  that 
propagation  is  in  the  x  direction,  with  y  being  the  horizontal  dimension  normal  to  propagation. 
The  domain  size  in y  that  has  been  used  is  300  wavelengths,  so  that  there  are  about  3.4  or  6.8 
grid  points  per  wavelength  of  sound.  (This  was  chosen  after  non-exhaustive  testing.)  The 
vertical  dimension  has  higher  resolution,  with  an  adjustable  aspect  ratio  A.  The  results  shown 
here  use  an  aspect  ratio  A  =  4.  Thus,  because  of  the  domain  mirroring  about  the  surface,  the 
vertical  domain  is  300/2/1  =  75  wavelengths.  The  source  can  be  positioned  at  any  depth  at  y  =  0. 


*  MATLAB'  is  a  high-level  language  and  interactive  environment  produced  by  The  Math  Works,  Inc. 
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For  the  results  shown  here  the  source  has  a  Gaussian  profile  in  z  andy  with  scale  length  of 
1 .25  m.  It  produces  an  axisymmetric  beam  with  a  width  of  about  20°.  More  recently  a  Green’s 
function  starting  field  has  been  implemented.  An  absorbing  boundary  is  required  to  prevent 
sound  wrap  around  and  aliasing  in  the  domain.  This  is  done  in  a  manner  similar  to  the  seafloor 
attenuation.  A  smooth  function  that  shrinks  from  unity  in  the  center  of  the  domain  to  near  zero 
at  the  edge  is  applied  at  each  edge  of  the  domain  as  a  factor  of  Ap.  The  function  for  each  edge  is 


where  8  is  the  step  length  in  the  propagation  (x)  direction,  Y  is  one/half  the  domain  width,  D  is 
the  distance  from  the  nearest  edge.  Parameter  values  used  successfully  with  400  Hz,  8  =  2k.  are 
a  -  kJ 50  and  /?=  2. 1 .  The  function  defined  over  the  entire  domain  Ap  is  composed  of  the  sum  of 
four  Aen  .  The  formula  is  written  to  be  compatible  with  the  Colosi  code  but  not  with  the 
Martin/Flatte  code,  which  implemented  a  radial  attenuation  function.  Successful  runs  without 
leakage  out  of  the  domain  (particularly  in  the  y  direction,  rather  than  z)  depend  on  careful  choice 
of  a  and  / 3 .  The  numbers  given  here  are  for  ay  =  300  k  domain,  400  Hz,  and  5  =  2k ,  and  may 
not  be  suitable  for  other  situations.  Algorithms  to  automatically  produce  suitable  parameters 
based  on  physical  principles  need  to  be  developed  and  implemented. 

Attenuation  in  the  seafloor  is  introduced  by  including  a  factor  Ag  in  the  medium-anomaly 
phase  fluctuations  at  each  step.  This  is  done  by  appropriately  reducing  the  magnitude  of  bottom- 
location  complex  phase-screen  elements  below  unity,  as  with  Ae.  That  is,  the  field  is  multiplied 
by  P-Ae  /fflexp(/(p)  at  each  step,  with  Ab  less  than  one  in  the  bottom  and  cp  determined  by  the 
fluctuating  sound-speed  of  the  medium. 

Propagation  is  modeled  with  domain  depths  of  60  m  and  80  m.  Background  sound  speed  in 
the  water  column  c(z)  has  a  three-layer  form,  with  fixed  dc/dz  in  each  layer.  The  sound  speed 
values  that  define  the  layers  are  c=1522  m/s  at  the  surface,  1520  at  15  m  depth,  1484  at  30  m, 
and  1481  at  the  bottom.  Seafloor  sound-speed  is  1650  m/s  at  the  interface,  increasing  with  a 
gradient  of  1  m/s  per  meter  below  the  seafloor.  Attenuation  in  the  seafloor  is  0.1  db/wavelength. 

Most  of  the  test  runs  include  time-frozen  nonlinear  wave-type  thermocline  displacements 
added  to  the  background  c(z).  These  have  the  form  of  permanent  wave  solutions  of  the 
Korteweg-de Vries  equation,  given  by  the  formula  r|  =  rjo  sech2(  (£-4o)  /  L  ) ,  where  E,  is  distance 
in  the  direction  of  travel  of  the  nonlinear  wave,  L  is  the  half  width,  and  rjo  is  the  maximum 
displacement  (the  amplitude).  As  appropriate  for  numerical  input,  wave  dimensions  are 
stretched  in  domain  slices  (across  the  wave)  that  lie  normal  to  the  propagation  direction. 

Run  time  is  1700  sec  for  a  range  of  10  km  (10242  FFT,  1333  steps,  9997.5  m)  on  a  2.8  GHz 
Pentium-4  Windows  PC  with  adequate  memory  (greater  than  1  GB).  Run  time  increases  if 
interim  fields  are  saved  or  more  complex  domains  are  input  or  computed  algorithmically  during 
the  run. 

Layered  Environment 

The  first  results  shown  here  will  be  for  a  uniform  layered  environment  with  no  internal 
waves.  Figures  1  and  2  show  the  amplitude  and  phase  of  sound  propagated  using  the  1024  x 
1024  point  option  in  the  60-m  domain.  Fields  at  ranges  of  500  and  5000  m  are  shown.  These 
show  the  effects  of  the  y-boundary  absorption  and  of  spreading.  The  results  appear  as  expected. 


3 


hr^JfdB) 


-500  -400  -300  -200  -100  0  100  200  300  400  500 

distance  (m) 

phase  angle 


500  400  300  200  100  0  100  200  300  400  500 

distance  (m) 


Figure  1.  The  acoustic  field  after  500-m  propagation  through  the  60-m  deep  three-layer 
environment  with  no  internal  waves  is  shown.  Intensity  in  the  tangential  cross-range  plane  is 
shown  at  the  top,  phase  in  the  same  plane  at  the  bottom.  The  beam  source  has  expanded  to  fill  a 
lateral  distance  of 200  m,  indicating  a  beam  width  of  about  22  degrees. 
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Figure  2.  Acoustic  field  results  for  5000-m  propagation  through  the  60-m  deep  three-layer 
environment.  This  range  is  approximately  twice  that  where  the  sound  impinges  on  the  lateral 

absorbing  layers. 
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Figure  3.  400-Hz  acoustic  field  after  propagation  of  20  km  through  the  80-m  deep  three-layer 
environment  in  2 A  steps  is  shown.  The  intensity  shows  a  ripple. 


Figure  4.  Intensity  and  phase  at  69.5  m  depth  from  the  previous  figure  are  plotted  as  functions 

of  cross-range  (y)  distance. 
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Figure  3  shows  results  for  20-km  propagation  in  the  deeper  domain.  At  this  range,  some 
spurious  focusing  occurs.  This  effect  can  be  caused  by  the  rather  large  range  step  size,  and  may 
also  emanate  from  the  edge  absorption  method.  Figure  4  is  a  line  drawing  of  amplitude  and 
phase  taken  from  Figure  3,  more  clearly  showing  the  amplitude  ripple. 

Internal  Wave  Environments 

Next,  Figure  5  shows  the  acoustic  field  after  propagation  through  an  environment  having  an 
internal  wave  aligned  to  be  consistent  with  motion  in  the  y  direction,  so  that  there  is  no  x 
dependence  of  the  sound  speed  field.  They’  dependence  of  c(x,y,z)  can  be  seen  in  the  figure. 

The  phase  plot  has  a  minimum  above  the  trough  of  the  wave,  indicating  a  maximum  speed,  as 
expected.  The  greater  intensity  to  the  right  of  the  wave  also  indicates  refraction  by  the  wave.  A 
meaningful  way  to  express  the  acoustic  effect  under  these  conditions  would  be  horizontal 
refraction  of  normal  modes,  with  group  speed  maxima  in  the  wave  expelling  the  mode  energy. 

The  fields  shown  thus  far  have  been  for  environments  with  no  x-direction  variability.  Now, 
results  are  shown  for  environments  with  range  variability  created  by  shifting  the  field  c(x=0  ,  y , 
z)  in  the  y-direction  a  fixed  number  of  grid  points  for  each  range  step,  creating  features  having 
fixed  angle  </>  with  respect  to  the  propagation  path  aty  =  0.  Figure  6  shows  results  for  a  field  of 
two  waves  that  move  across  and  out  the  domain  in  the  y  direction,  then  re-entering  the  domain. 
The  wave  field  rj(x,y)  is  shown  in  Figure  7. 
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Figure  5.  The  acoustic  field  at  distance  5000  m  from  the  source  is  shown  for  the  situation  of 
propagation  parallel  to  the  crest  of  an  internal  wave  with  L=86  m  and  rjo  =  12  m.  The  60-m 
depth  environment  is  used.  The  boundaries  of  the  two  layers  are  drawn  in  this  view  looking  in 
the  direction  of propagation  (x),  showing  the  fixed  position  of  the  internal  wave  in  the  cross - 
range  direction  (y).  The  source  is  at  y  -  0  and  depth  40  m.  Note  the  minimum  of  phase  above 
the  trough  of  the  wave,  and  the  larger  intensity  on  the  right  caused  by  refraction  away  from  the 

warm  anomaly  in  the  wave. 
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Figure  6.  The  field  resulting  from  propagation  through  internal  waves  of  10-m  amplitude  and 
95-m  half  width  (L)  at  an  angle  of  <f>  =  23.7  °.  This  is  oblique  acoustic  propagation  with  respect 
to  the  wave  troughs.  Because  the  down-range  half-width  of  the  internal  waves  is  235  m,  there  is 
a  strong  possibility  that  coupled  mode  propagation  is  occurring.  Adiabatic  mode  refraction  of 
the  type  illustrated  in  Figure  5  is  not  expected  to  occur. 


y(m) 


Figure  7.  The  internal  wave  environment  used  to  generate  the  field  of  Figure  6  is  shown.  The 
depth  of  the  top  of  the  strong  gradient  layer  is  plotted  versus  x  and  y.  Wave  troughs  are  blue. 
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Figures  8-11  show  results  for  simulations  having  one  internal  wave  of  r|0=10  m,  L= 95  m  at 
4950  m  from  the  source  at  >>  =  0,  with  the  wave  at  various  angles  <j>.  These  are  similar  to  the 
situation  shown  in  Figure  7  except  only  one  wave  is  present.  The  fields  that  result  at  the  range 
of  x  =  10  km  (Figures  8  and  9)  have  significant  structure  in  y,  as  expected  from  either  coupled 
mode  or  adiabatic  propagation  through  media  that  have  3-D  variability. 

Figure  10  shows  the  depth-integrated  range -normalized  intensity  (intensity  times  range)  that 
results  from  waves  at  each  of  the  two  small  angles  of  <j>  equal  to  4°  and  8°.  The  smaller  angle 
case  clearly  shows  refraction,  the  other  does  not,  although  the  next  figure  will  suggest  that 
refraction  also  occurs.  Figure  1 1  shows  the  differences  between  the  fields  of  Figure  10  and 
control  cases  with  no  internal  waves.  In  the  top  panel,  the  blue  area  below  the  4°  internal-wave 
position  has  depleted  energy,  with  energy  diverted  into  the  beam  near  the  wave  seen  at  the  top. 
The  same  effect  can  be  seen  in  the  lower  panel  for  the  8°  wave.  Each  of  the  panels  also  shows 
azimuthal  energy  dependence  appearing  as  a  few  beams.  The  beams  are  green,  separated  by  blue 
areas  of  energy  depletion.  This  effect  is  caused  by  source  to  internal-wave  distance  being  a 
function  of  azimuth,  causing  azimuthal  variations  in  mode  coupling,  and  subsequent  variations  in 
energy  loss  via  bottom  interaction  [7]. 
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Figure  8.  The  intensity  and  phase  at  the  10-ktn  range  end-plane  for  the  case  of  one  internal 
wave  5  km  from  the  source  (measured  at  y  -  0).  The  internal  wave  is  at  angle  <j>  =  16° .  This  is 
a  situation  that  shows  no  evidence  of  horizontal  deflection.  The  phase  minimum  is  at  the  point 
closest  to  the  source,  which  is  y  =  0  in  the  center  of  the  figure. 
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Figure  9.  Intensity  and  phase  are  shown  for  a  situation  similar  to  that  of  the  previous  figure, 
except  the  wave  angle  is  <f>=  4°.  This  situation  shows  horizontal  refraction  of  the  sound.  The 
depth-averaged  sound  intensity  is  not  symmetric  around  y  —  0,  unlike  the  <f>—16  °case  (not 

shown).  The  phase  minimum  is  not  at  y  =  0. 


Figure  10.  The  depth-averaged  range-normalized  intensity  is  shown  in  plan  view  for  two  cases, 
each  with  a  single  internal  wave  4950  m  from  the  source  at  y  =  0.  (Upper  panel)  The  wave  has 
</>  =  4  °.  (Lower  panel)  The  wave  has  <f>  =  8°.  The  trough  of  the  wave  is  shown  in  each  panel  in 
white.  The  effects  of  refraction  are  clearly  seen  in  the  upper  panel.  [The  8  °  wave  situation 
shows  slight  azimuthal  dependence  of  the  energy,  the  result  of  weak  horizontal  refraction  and 
possibly  also  altered  bottom  attenuation  effects  resulting  from  mode  coupling.] 
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Figure  11.  The  differences  of  the  fields  of  the  previous  figure  with  the  no  internal  wave  case  are 
plotted.  The  intensity  scales  differ  in  the  two  plots.  Amplification  along  the  wave  caused  by 
refraction  is  seen  in  both  cases.  Positive  values  in  the  refracted  beam  in  the  upper  image  are 

clipped. 


Figure  12.  Two  related  environments  are  shown  in  the  format  of  Figure  7.  These  are  identical 
along  the  line  y  =  0.  On  the  left  is  a  strongly  azimuthally  varying  field  with  a  single  wave  at 
angle  </>  -8°;  on  the  right  is  the  weakly  azimuthally  varying  case  generated  by  using  they  -  0 

field  at  all  y. 
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angle:  41 .3263’  angle:  36  2332’ 


angle  23.7335° 


Figure  13.  3-D  to  pseudo-2-D  comparison  for  various  angles  of  waves.  At  the  final  range  of x  - 
10  km  and y  -  0,  intensity  is  shown  above,  phase  below.  The  internal  waves  all  have  amplitude 
10  m,  horizontal  scale  L  =  95  m.  The  most  pronounced  discrepancy  that  is  seen,  at  <f>  =23°  is 

unexplained. 


angle:  16.3365* 


-2  0  2 
Phase  (radians) 


angle:  8.3377* 


Phase  (radians) 


angle:  4.19* 


-2  0  2 
Phase  (radians) 


Figure  14.  Additional  3-D  to  pseudo-2-D  comparisons  for  10:95  waves.  The  <f>  =  4  °case  shows 
a  clear  difference,  indicating  refractive  effects,  the  others  do  not. 
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3-D  to  2-D  comparison 

The  accuracy  of  A^x  2D  simulations  of  propagation  through  3D  environments,  such  as  in 
previous  work  [7],  can  be  evaluated  using  the  3D  model.  To  avoid  discrepancies  cause  by 
comparing  codes,  we  perform  such  an  evaluation  by  comparing  results  from  strongly-varying 
and  weakly  azimuthally  varying  3D  environments  that  are  identical  along  the  liney  =  0.  Two 
environments  of  this  type  are  shown  in  Figure  12. 

Using  environment  pairs  of  the  type  shown  in  Figure  12  the  3D  to  2D  comparison  is  made 
for  single  10-m  amplitude,  L  =  95  m  waves  of  depression  propagating  at  various  angles.  Figures 
13  and  14  show  that  the  2D  simulations  give  very  accurate  results  for  internal  waves  having 
crests  (or  troughs)  at  angle  8  degrees  or  greater  from  the  direction  of  propagation.  These  are  the 
same  situations  that  show  little  refractive  influence  on  intensity. 

Mode  Coupling 

The  process  of  coupled  mode  propagation  has  been  demonstrated  to  occur  when  sound  in 
waveguides  similar  to  that  modeled  here  encounters  internal  waves  of  the  type  used  here.  Figure 
15  shows  that  this  is  indeed  occurring  for  a  wave  at  angle  tf>=  30.4°.  The  modeled  wave  has  a  = 

1 0  m  and  L  =  95  m  as  previously,  which  has  a  wave  scale  of  1 88  m  along  y  =  0.  This  case  is 
expected  to  have  mode  coupling  [8],  verified  by  the  modal  decomposition  of  the  field  as  a 
function  of  x  shown  in  the  figure. 


1000  2000  3000  4000  5000  6000  7000  8000 

range  (m) 

Figure  15.  The  magnitudes  of  mode  coefficients  times  range1  2  are  computed  along  y  =  0 
for  a  case  similar  to  those  shown  earlier  (400  Hz,  source  depth  40  m,  one  10:95  nonlinear 
internal  wave).  The  internal  wave  is  at  angle  30.4  °with  position  x  =  4000  m  aty  =0.  Mode 
filtering  is  performed  with  mode  functions  (computed  using  the  code  FEMODE.F  from  Naval 
Research  Lab.)  for  the  unperturbed  environment,  i.e.  outside  the  wave.  The  area  within  the  wave 
where  these  are  not  the  correct  mode  functions  is  shown  in  gray.  Modes  should  be  uncoupled 
outside  the  wave,  so  the  oscillations  in  range  of  the  coefficients  show  small  numerical  errors 
known  to  be  associated  with  overly  long  range  step  increment.  Despite  this,  mode  coupling  in 
the  wave,  which  is  expected,  can  be  seen  to  strip  energy  from  modes  5  and  6. 
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Figure  16.  A  result  similar  to  that  of  the  previous  figure  is  shown,  with  the  exception  that  the 
wave  angle  is  16  degrees.  The  internal  wave  scale  projected  onto  a  radial  line  of  propagation  at 

y  =  0  is  Lp-  95/sin<j>  —  344  m. 


For  decreasing  wave  angle  ^the  scale  of  the  internal  waves  projected  onto  radial  lines  of 
propagation  becomes  quite  large,  which  reduces  lateral  gradients  and  weakens  the  mode 
coupling  effect  [8],  Figure  16  shows  that  despite  this  trend,  a  wave  at  angle  16  degrees  which 
has  an  effective  wave  scale  near  350  m  along  the  propagation  path  still  causes  significant 
coupling,  with  mode  5  severely  depleted  after  traversal  of  the  wave. 

Vertical  Channel:  PIES  surface  reflection  simulation 

Propagation  of  high-frequency  sound  (6  to  15  kHz)  emitted  from  a  bottom-mounted  PIES 
instrument  to  the  ocean  surface  and  then  back  to  the  bottom  has  also  been  modeled.  A  uniform 
sound  speed  has  been  used  throughout,  enabling  an  analytic  field  to  be  applied  at  the  surface, 
reducing  the  computation  by  50  percent.  Upward-propagation  test  runs  reproduced  the  phase 
calculated  using  spherical  spreading  with  high  accuracy.  The  surface  reflection  is  done  with  a 
phase-screen  treatment;  at  each  point,  the  field  is  multiplied  by  P  =  exp(/cp),  with  phase  of  (p  = 

4 nhtA.  The  test  simulations  typically  used  32 -X  steps,  400-X  domain  widths  for  the  1024  x  1024 
point  grid,  and  800-A.  domains  for  2048  x  2048  point  grids.  The  edge  attenuation  functions  have 
different  coeffecients  from  the  low-frequency  case,  with  P=8/8  for  the  wide  case  with  more  grid 
points,  double  its  value  for  the  other  case.  The  value  a  =  X/2  is  used. 

Figure  17  shows  the  surface  wave  field  that  was  used.  This  is  computed  by  starting  with  a 
Pierson-Moskowitz  spectrum  with  noise  added,  also  having  an  azimuthal  taper.  This  is  Fourier 
transformed  into  the  spatial  domain.  The  field  is  self-similar  at  the  wavelengths  near  the  Fresnel 
scale  that  are  important  to  the  acoustics  ( R/=  DA  where  D  is  depth),  so  the  field  is  merely 
stretched  to  fit  domains  of  different  sizes  for  these  test  runs. 
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Figure  1 7.  The  gravity  wave  field  at  the  surface  used  in  the  PIES  vertical  channel 

simulations. 


Fields  that  result  at  the  bottom  in  512-m  depth  for  12  kHz  are  shown  in  Figure  18.  Two 
different  step  sizes  are  used,  showing  restrictions  on  the  SSF  PE  in  a  narrow  domain.  Short  steps 
are  required  to  avoid  aliasing.  This  is  not  a  problem  for  the  shallow-water  waveguide  because 
short  steps  (relative  to  wavelength)  are  needed  to  model  the  highly-variable  medium  in  the  lower 
frequency  case,  shorter  than  would  otherwise  be  needed.  The  vertical  channel  can  be  relatively 
uniform,  but  this  can  not  be  taken  advantage  of  with  very  long  steps  because  of  the  aliasing. 
Considering  the  acoustic  field,  note  the  intense  focus  (with  interference  pattern)  in  the  lower  part 
of  Figure  18  that  appears  beneath  the  concave  (from  below)  surface  wave  crest.  Figures  19  and 
20  show  interesting  patterns  at  10  kHz  for  water  depths  of  3000  and  512  m,  respectively. 

Future  Work 

A  number  of  enhancements  would  improve  the  code  and  would  make  simulation  results 
more  suitable  for  research  purposes.  These  are  listed  here. 

1 .  The  density  must  be  allowed  to  change.  This  can  be  done  by  changing  the  operator  P 
which  is  applied  in  the  spatial  domain.  This  has  been  implemented  but  not  fully  tested. 
Methods  for  this  are  discussed  in  the  literature  [9], 

2.  Theoretical  restrictions  for  the  grid  spacing  and  the  step  length  should  be  computed  as  a 
function  of  frequency,  and  user-generated  grids  should  be  tested  for  compliance. 

3.  An  algorithm  is  needed  to  automatically  generate  an  effective  edge  absorption  function, 
defined  in  terms  of  the  grid  spacing  and  the  step  increment.  An  entirely  new  function 
might  be  in  order. 

With  these  improvements  the  code  will  give  more  reliable  results  and  will  be  suitable  for 
studying  3D  propagation  in  arbitrarily  complicated  environments.  The  code  can  already  handle 
any  water  column  and  bathymetric  conditions,  although  the  input  scheme  for  complicated 
environments  has  not  been  settled  upon. 
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Figure  18.  Edge  absorption  test  for  the  vertical  propagation  channel.  Results  are  compared  for 
128  4-m  steps  (top)  versus  4  128-m  steps  (bottom).  Amplitude  is  shown  on  the  left,  phase  on  the 
right.  The  propagator  moves  sound  through  the  boundaries  in  the  long-step  case,  causing 

aliasing. 
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Figure  19.  Acoustic  amplitude  (top)  and  phase  (bottom)  for  a  10-kHz  simulation  of  a  surface 
reflected  field from  a  bottom-mounted  PIES  in  3000-m  depth  water.  2048  jc  2048  domain. 
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10  kHz.  Wide  angle  PE.  2048  2 . 512  m  depth 


Figure  20.  Intensity  resulting  from  a  wide-angle  PE  run  vertically  upward  and  back  in  512  m 
depth  water.  The  frequency  is  10  kHz,  and  the  surface  waves  are  as  shown  in  Figure  1 7. 
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